Enzyminator 3000

---

Overview

All proccesses in cells rely on chemical reactions, which can require lots of energy. The Second Law of Thermodynamics suggests that entropy, or disorder, naturally tends to increase, but life itself seems to contradict that, with it definiteionallyy being the organization of matter in a self-sustaining system.

Here's where enzymes come it. Enzymes are vital catalysts in these cells to lower the activation energy needed for biological processes. As such, when producing a large quantity of a product or needing to produce a response quickly, enzymes are almost always used to . The program modularly mimics the chaotic "random" collision in a cell to model how long it takes for a substrate to reach an enzyme.

---

The Program

The following properties will be variable:

• Number of particles
• Size (radius) of particles
• Mass of particles
• Shape of Boundary
• Size of Boundary
• Temperature (Avg. Kinetic Energy)

In the first mode, these particles will be added in a random configuration of positions and velocity directions, with an animation plyaing out until an enzyme and substrate collide and a reaction occurs. The number of steps will then be displayed.

In the second mode, the program goes through multiple trials with varying numbers of substrate particles. The program averages these out and plots the number of substrate particles against the reaction time in a bar graph with error bars. It should theoretically show an inverse relationship as the reaction time decreases as substrate particles increase.

import sys
import math
import random as random
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import *
from matplotlib.patches import *  # Used to draw rectangle
random.seed(None)        # Seed generator, None => system clock

stat = input("Stat mode (y/n)? ")
if stat=='y':
    stat = True
elif stat=='n':
    stat = False
else:
    stat = False

if stat:
    %matplotlib inline
else:
    %matplotlib notebook

class Particle:
    def __init__(self):
        self.radius = 0
        self.mass = 0
        self.velocity = 0
        self.x = 0
        self.y = 0
        self.angle = 0
        
def colliding(p1, p2):
    return math.hypot(p1.x-p2.x, p1.y-p2.y) <= (p1.radius + p2.radius)

def colliding_n(p1, p2):
    return math.hypot(p1.nx-p2.nx, p1.ny-p2.ny) <= (p1.radius + p2.radius)

class RectangleBoundary:
    name = "rectangle"
    property_keys = ["width", "height"]
    
    def __init__(self):
        self.properties = dict.fromkeys(self.property_keys)

    def draw(self):
        plt.gca().add_patch(Rectangle((-self.properties["width"]/2,-self.properties["height"]/2), self.properties["width"], self.properties["height"], fill=False, edgecolor='k',lw=2))
        plt.xlim([-self.properties["width"]/2-3,self.properties["width"]/2+3])
        plt.ylim([-self.properties["height"]/2-3,self.properties["height"]/2+3])
    
    def check_inside(self, x, y, r):
        return (-self.properties["width"]/2+r < x < self.properties["width"]/2-r) and (-self.properties["height"]/2+r < y < self.properties["height"]/2-r)

    def get_angle(self, x, y, r):
        if(y>=self.properties["height"]/2-r and x>=self.properties["width"]/2-r):
            return 3*np.pi/4
        elif(y<=-self.properties["height"]/2+r and x<=-self.properties["width"]/2+r):
            return 3*np.pi/4
        elif(y<=-self.properties["height"]/2+r and x>=self.properties["width"]/2-r):
            return 5*np.pi/4
        elif(y>=self.properties["height"]/2-r and x<=-self.properties["width"]/2+r):
            return 5*np.pi/4
        elif (y>=x*(self.properties["height"]-2*r)/(self.properties["width"]-2*r)) == (y>=-x*(self.properties["height"]-2*r)/(self.properties["width"]-2*r)):
            return 0
        else:
            return np.pi/2
    
    def get_limits(self):
        return ((-self.properties["width"]/2, self.properties["width"]/2), (-self.properties["height"]/2, self.properties["height"]/2))

class CircleBoundary:
    name = "circle"
    property_keys = ["diameter"]
    
    def __init__(self):
        self.properties = dict.fromkeys(self.property_keys)

    def draw(self):
        plt.gca().add_patch(Circle((0, 0), self.properties["diameter"]/2, fill=False, edgecolor='k',lw=2))
        plt.xlim([-self.properties["diameter"]/2-3,self.properties["diameter"]/2+3])
        plt.ylim([-self.properties["diameter"]/2-3,self.properties["diameter"]/2+3])
    
    def check_inside(self, x, y, r):
        return x**2 + y**2 < (self.properties["diameter"]/2 - r)**2

    def get_angle(self, x, y, r):
        return np.arctan2(-x,y)%(2*np.pi)
    
    def get_limits(self):
        return ((-self.properties["diameter"]/2, self.properties["diameter"]/2), (-self.properties["diameter"]/2, self.properties["diameter"]/2))
    
boundaries = [RectangleBoundary(), CircleBoundary()]
boundaries_string = ", ".join([b.name for b in boundaries])
boundary = 0
while boundary==0:
    response = input(f"Select a boundary type ({boundaries_string}): ")
    for b in boundaries:
        if b.name == response:
            boundary = b
            break
for key in boundary.property_keys:
    boundary.properties[key] = float(input(f"Enter {key}: "))
limits = boundary.get_limits()
print()

if(stat):
    substrates = []
else:
    substrates = [Particle() for i in range(int(input("Enter number of substrates: ")))]
sradius = float(input("Enter substrate radius: "))
for s in substrates:
    s.radius = sradius
smass = float(input("Enter substrate mass: "))
for s in substrates:
    s.mass = smass
for i in range(len(substrates)):
    s = substrates[i]
    while(True):
        s.x = random.uniform(limits[0][0], limits[0][1])
        s.y = random.uniform(limits[1][0], limits[1][1])
        if(not boundary.check_inside(s.x, s.y, s.radius)):
            continue
        repeat = False
        for j in range(i):
            if(colliding(s, substrates[j])):
                repeat = True
                break
        if(repeat):
            continue
        break
print()

enzymes = [Particle() for i in range(int(input("Enter number of enzymes: ")))]
radius = float(input("Enter enzyme radius: "))
for e in enzymes:
    e.radius = radius
mass = float(input("Enter enzyme mass: "))
for e in enzymes:
    e.mass = mass
for i in range(len(enzymes)):
    e = enzymes[i]
    while(True):
        e.x = random.uniform(limits[0][0], limits[0][1])
        e.y = random.uniform(limits[1][0], limits[1][1])
        if(not boundary.check_inside(e.x, e.y, e.radius)):
            continue
        repeat = False
        for j in range(i):
            if(colliding(e, enzymes[j])):
                repeat = True
                break
        for s in substrates:
            if(colliding(e, s)):
                repeat = True
                break
        if(repeat):
            continue
        break
print()

solvents = [Particle() for i in range(int(input("Enter number of solvents: ")))]
radius = float(input("Enter solvent radius: "))
for s in solvents:
    s.radius = radius
mass = float(input("Enter solvent mass: "))
for s in solvents:
    s.mass = mass
for i in range(len(solvents)):
    s = solvents[i]
    while(True):
        s.x = random.uniform(limits[0][0], limits[0][1])
        s.y = random.uniform(limits[1][0], limits[1][1])
        if(not boundary.check_inside(s.x, s.y, s.radius)):
            continue
        repeat = False
        for j in range(i):
            if(colliding(s, solvents[j])):
                repeat = True
                break
        for su in substrates:
            if(colliding(s, su)):
                repeat = True
                break
        for e in enzymes:
            if(colliding(e, s)):
                repeat = True
                break
        if(repeat):
            continue
        break
print()

KE = float(input("Enter temperature (average kinetic energy prop. to speed^2): "))

def setup_points():
    particles = substrates+enzymes+solvents
    for i in range(len(particles)):
        p = particles[i]
        while(True):
            p.x = random.uniform(limits[0][0], limits[0][1])
            p.y = random.uniform(limits[1][0], limits[1][1])
            if(not boundary.check_inside(p.x, p.y, p.radius)):
                continue
            repeat = False
            for j in range(i):
                if(colliding(p, particles[j])):
                    repeat = True
                    break
            if(repeat):
                continue
            break
    for p in particles:
        p.velocity = np.sqrt(KE/p.mass)
        p.angle = random.uniform(0, 2*np.pi)


circles = []

def init():
    for s in substrates:
        s.circle = ax.add_patch(plt.Circle((s.x, s.y), s.radius, color='r'))
        circles.append(s.circle)
    for e in enzymes:
        e.circle = ax.add_patch(plt.Circle((e.x, e.y), e.radius, color='g'))
        circles.append(e.circle)
    for s in solvents:
        s.circle = ax.add_patch(plt.Circle((s.x, s.y), s.radius, color='b'))
        circles.append(s.circle)
    return circles

def step(frame):
    for i in range(len(substrates+enzymes+solvents)):
        p = (substrates+enzymes+solvents)[i]
        
        p.nx = p.x + p.velocity*np.cos(p.angle)
        p.ny = p.y + p.velocity*np.sin(p.angle)

        if(not boundary.check_inside(p.nx, p.ny, p.radius)) :
            surface_angle = boundary.get_angle(p.nx, p.ny, p.radius)
            p.angle = 2*surface_angle - p.angle
            p.nx = p.x
            p.ny = p.y
            counter=0
            while(counter < 10 and not boundary.check_inside(p.nx, p.ny, p.radius)):
                p.nx += p.velocity*np.cos(p.angle)
                p.ny += p.velocity*np.sin(p.angle)
                counter+=1
            if counter>=10:
                p.nx = p.x
                p.ny = p.y
        
        for q in (substrates+enzymes+solvents)[:i]:
            if colliding_n(p, q):
                if((p in substrates and q in enzymes) or (q in substrates and p in enzymes)):
                    if stat:
                        return True
                    pprop = p.radius/(p.radius+q.radius)
                    qprop = q.radius/(p.radius+q.radius)
                    ax.plot(pprop*q.x+qprop*p.x, pprop*q.y+qprop*p.y, 'y*', ms=15)
                    ax.text(0.5,0.5,'# of steps =' + str(frame), ha='center', va='center',fontsize=20, transform=ax.transAxes)
                    anim.event_source.stop()
                c = np.arctan2(p.ny-q.ny, p.nx-q.nx)
                vpx = np.cos(c)*(p.velocity*np.cos(p.angle-c)*(p.mass-q.mass)+2*q.mass*q.velocity*np.cos(q.angle-c))/(p.mass+q.mass) + p.velocity*np.sin(p.angle-c)*np.cos(c+np.pi/2)
                vpy = np.sin(c)*(p.velocity*np.cos(p.angle-c)*(p.mass-q.mass)+2*q.mass*q.velocity*np.cos(q.angle-c))/(p.mass+q.mass) + p.velocity*np.sin(p.angle-c)*np.sin(c+np.pi/2)
                vqx = np.cos(c)*(q.velocity*np.cos(q.angle-c)*(q.mass-p.mass)+2*p.mass*p.velocity*np.cos(p.angle-c))/(q.mass+p.mass) + q.velocity*np.sin(q.angle-c)*np.cos(c+np.pi/2)
                vqy = np.sin(c)*(q.velocity*np.cos(q.angle-c)*(q.mass-p.mass)+2*p.mass*p.velocity*np.cos(p.angle-c))/(q.mass+p.mass) + q.velocity*np.sin(q.angle-c)*np.sin(c+np.pi/2)
                p.velocity = math.hypot(vpx,vpy)
                p.angle = np.arctan2(vpy,vpx)
                q.velocity = math.hypot(vqx,vqy)
                q.angle = np.arctan2(vqy,vqx)
                p.nx = p.x
                p.ny = p.y
                counter = 0
                while(counter < 10 and colliding_n(p, q)):
                    p.nx += p.velocity*np.cos(p.angle)
                    p.ny += p.velocity*np.sin(p.angle)
                    counter+=1
                if counter>=10:
                    p.nx = p.x
                    p.ny = p.y
                
    for p in substrates+enzymes+solvents:
        p.x = p.nx
        p.y = p.ny
        if not stat:
            p.circle.center = (p.x, p.y)
            
    if stat:
        return False
    
    return circles

fig, ax = plt.subplots()

if(stat):
    x = []
    y = []
    error = []
    
    for i in range(1,15):
        x.append(i)
        p = Particle()
        p.mass = smass
        p.radius = sradius
        substrates.append(p)
        
        t = []
        
        for j in range(200):
            setup_points()
            frame = 0
            while True:
                frame+=1
                s = step(frame)
                if(s):
                    t.append(frame)
                    break
        y.append(np.mean(np.array(t)))
        error.append(2*np.std(np.array(t))/np.sqrt(len(t)))
    
    ax = fig.add_axes([0,0,1,1])
    ax.bar(x, y, yerr=error)
    
    plt.title('Enzyminator')
    plt.xlabel("Number of Substrates")
    plt.ylabel("Average Number of Steps")
    
    plt.show()
            
else:
    setup_points()
    
    ax.set_xlim([limits[0][0]-3,limits[0][1]+3])
    ax.set_ylim([limits[1][0]-3,limits[1][1]+3])
    boundary.draw()

    plt.title('Enzyminator')
    plt.axis("equal")

    Writer = writers['ffmpeg']
    writer = Writer(fps=15, metadata=dict(artist='Me'), bitrate=1800)

    anim = FuncAnimation(fig, step, init_func=init, interval=20, blit=True, repeat=False)

---

Conclusion

Here is a demonstration of the Brownian motion in a cell:

With all the water molecules in the way, it is very inefficient, taking lots of time to react. If we increase the temperature by 4 times (effectively doubling the initial velocity), the reaction rate does increase.

However, increasing temperature is not a good idea in cells because these enzymes only work in a narrow range of temperatures and will denature otherwise. Instead, cells just produce many redundant copies of substrates or enzymes to increase the reaction rate. Here, we've added 2 more substrates.

To investigate the effect of the number of substrates of reaction rate, we varied the number of substrate particles while keeping all else constant and measured the average number of steps for one of them to meet an enzyme. Here is the data collected:

The graph shows a clear inversely trend. As the number of substrate particles increases demonstrates an increase in reaction rate, the average reaction time decreases. The time decreases rapidly at first but eventually tapers off, as demonstrated by the 95% error bars, which are distinct between trials at lower numbers but overlap at higher levels. This makes sense as an increase from 1 to 2 is a much larger relative change when compared to 13 to 14. In addition, adding more and more substrates leads to diminishing returns as the reaction rate's limiting factor becomes the number of available enzymes instead of substrates.

To conclude, the model shows that as the number of substrates increases, so does the reaction rate of enzymes.